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("^ , EPW (Electron-Phonon coupling using Wannier functions) is a program written in F0RTRAN90 for 

^vj ■ calculating the electron-phonon coupling in periodic systems using density-functional perturbation 

theory and maximally-localized Wannier functions. EPW can calculate electron-phonon interaction 
self-energies, electron-phonon spectral functions, and total as well as mode-resolved electron-phonon 
coupling strengths. The calculation of the electron-phonon coupling requires a very accurate sam- 
pling of electron-phonon scattering processes throughout the Brillouin zone, hence reliable calcula- 
■ tions can be prohibitively time-consuming. EPW combines the Kohn-Sham electronic eigenstates and 

^^ , the vibrational eigenmodes provided by the Quantum-ESPRESSO package [ij with the maximally local- 

ized Wannier functions provided by the waimier90 package 2] in order to generate electron-phonon 
matrix elements on arbitrarily dense Brillouin zone grids using a generalized Fourier interpolation. 
5— ( , This feature of EPW leads to fast and accurate calculations of the electron-phonon coupling, and 

^~^ • enables the study of the electron-phonon coupling in large and complex systems. 
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PACS numbers: 63.20.kd, 63.20.kk, 71.15.-m, 74.25.Kc, 74.70.-b 



PROGRAM SUMMARY 

^ : Program Title: EPW 
■ O Journal Reference: 

R , Catalogue identifier: 

Licensing provisions: GNU Public License 
Programming language: Fortran 90 

Computer: any architecture with a Fortran 90 compiler 
►^ RAM: heavily system dependent, as small as a few MB 

ryr\ Number of processors used: optimized for 1 to 64 processors 

7— H Classification: 7 

■^ External routines /libraries: BLAS, LAPACK, MPI, FFTW 

\l ,, Subprograms used: wannier90 

l/~) ' Nature of problem: The calculation of the electron-phonon coupling from first principles requires a very accurate sampling 
^^ ,, of electron-phonon scattering processes throughout the Brillouin zone; hence reliable calculations can be prohibitively time- 
consuming. 

Solution method: EPW makes use of a real-space formulation and combines the Kohn-Sham electronic eigenstates and the 
vibrational eigenmodes provided by the Quantum-ESPRESSO package with the maximally localized Wannier functions provided 
by the wannier90 package in order to generate electron-phonon matrix elements on arbitrarily dense Brillouin zone grids using 
» 1 ' a generalized Fourier interpolation. 
Cu ^ Running time: single processor examples typically take 5-10 minutes 

References: Giustino, F., Cohen, M. L., and Louie, S. G., Physical Review B 76 (2007) 165108. 

I. INTRODUCTION 

The electron-phonon interaction plays a crucial role in the electron and lattice dynamics of condensed matter 
systems. For example, phenomena such as the electrical resistivity [3| and conventional superconductivity [J] are 
a direct consequence of the interaction between electrons and lattice vibrations. The electron-phonon interaction 
also plays an important role in the thermoelectric effect J5|. Other fundamental physical phenomena such as the 
Kohn effect [6| and the Peierls [7| distortions are also direct consequences of the electron-phonon interaction. The 
electron-phonon interaction is also responsible for the broadening of the spectral lines in angle-resolved photoemission 
spectroscopy [8| and in vibrational spectroscopies |9| , as well as for the temperature dependence of the band gaps in 



semiconductors [lOJ. 

The calculation of the electron-phonon coupling from first-principles is challenging because of the necessity of 
evaluating Brillouin zone integrals with high accuracy. Such calculation requires the evaluation of matrix elements 
between electronic states connected by phonon wavevectors [11[. Well-established software packages are available for 
computing electronic states and eigenvalues through density- functional theory [l|,ll^ll3|, as well as phonon frequencies 
and eigenmodes through density- functional perturbation theory IJ] . However large numbers of matrix elements may 
be necessary to achieve numerical convergence of the Brillouin zone intergrals over these matrix elements. 

For example, in order to compute the electronic lifetimes associated with the electron-phonon interaction it is 
necessary to evaluate a Brillouin zone integral over all the possible phonon wavevectors (thousands to millions). Since 
lattice-dynamical calculations for each phonon wavevector are at least as expensive as self-consistent total energy 
minimizations, achieving numerical convergence in the Brillouin zone integrals over the phonon wavevectors by brute- 
force calculations may become a prohibitive computational task. EPW exploits the real-space localization of Wannier 
functions to generate large numbers of first-principles electron-phonon matrix elements through a generalized Fourier 
interpolation. EPW therefore enables affordable, accurate, and extremely efhcient calculations of the electron-phonon 



coupling 15]. The use of maximally localized Wannier functions fMLWFs ) |16l .ll7| to calculate Brillouin zone integrals 
with high accuracy has been the object of a number of other studies |18l430l |. 

In this communication we outline the functionalities of EPW and the details of the technical release (Sec. II), we 
review the individual sections of the code (Sec. Ill), and we describe its parallel structure (Sec. VI). We then illustrate 
the capabilities of EPW through example calculations (Sec. V). We finally discuss some future directions for additional 
development (Sec. VI). Some technical details are given in the Appendix (Sec. VII). 

II. FUNCTIONALITIES AND TECHNICAL RELEASE 

EPW is a program written in F0RTRAN90 which calculates the electron-phonon coupling from first principles using 
maximally localized Wannier functions. EPW uses information provided by Qusintum-ESPRESSO [Ij and wannier90 [2], 



and runs as a post-processing tool. Electrons are described using density functional theory (DFT) [12|, |13| with plane- 
waves and pseudopotentials, cither separable norm-conserving |3ll - l34| or Vanderbilt ultrasoft [35!|. Lattice dynamical 
properties are calculated within density functional perturbation theory (DFP T) [14] . The theoretical background and 
methodology are thoroughly described in Ref. [l5|. Examples of quantities [36| which can be computed using EPW 
include: 

• the total electron-phonon coupling strength A, 

• the phonon self-energy associated with the electron-phonon interaction within the Migdal approximation, 

• the electron self-energy associated with the electron-phonon interaction within the Migdal approximation, 

• the Eliashberg electron-phonon spectral function a^F , 



• 



the transport electron-phonon spectral function a F^ 



EPW is freely available under the GNU General Public License (GPL). The current version is developed and main- 
tained using Subversion and is accessible to prospective developers and end-users at the website epw.org.uk. EPW 
employs the freely available FFTW, BLAS, LAPACK libraries in conjunction with several subroutines distributed within 
the Quantum-ESPRESSO package. Several subroutines in EPW are based upon modified Quantum-ESPRESSO subroutines 
as permitted under the GPL. The parallelization is achieved through the MPI library specification for message passing. 
The current version of EPW, v2.3, includes approximately 9000 independent lines of FQRTRAN90 code. In addition to 
the source code, several complete example calculations are provided with the EPW distribution. 

The inputs to EPW are as follows: 

• Phonon dynamical matrices for the wavevectors of a uniform Brillouin-zone grid centered at F (prefix . dyn files 
from Quantum-ESPRESSO). Only wavevectors within the irreducible wedge of the Brillouin zone are required. 

• The derivatives of the self-consistent potential with respect to the phonon perturbations, for the same wavevec- 
tors as above (prefix. dvscf files from Quantum-ESPRESSO). 

• The electron eigenfunctions and eigenvalues for the wavevectors of a uniform Brillouin-zone grid centered at F 
(prefix. wfc or prefix. save/ files from Quantum-ESPRESSO). 



• Norm-conserving pseudopotentials 32^ or Vanderbilt ultrasoft pseudopotentials 35 [. 

• A plain text input file specifying the runtime parameters. 

The electronic wavefunctions are calculated on a uniform grid using Quantum-ESPRESSD. Dynamical matrices and the 
derivatives of the self-consistent potential are also computed within Quantum-ESPRESSO for phonons in the irreducible 
wedge of the Brillouin zone. When choosing initial electron and phonon grids, it is necessary that the Brillouin 
zone grid for phonons be comensurate with the Brillouin zone grid for electrons in order to map the wavefunctions 
'0mk+q(r) onto V'mk'+Glr), with G a reciprocal lattice vector. For example, if the calculation is performed using a 
Brillouin zone grid of size 6x6x6 for the phonons, then the natural choices for the electronic Brillouin zone grid 
are either 6x6x6 or 12x12x12. This does not represent a computational bottleneck as phonon calculations are 
considerably more time-consuming than non-selfconsistent electronic calculations. 

III. COMPUTATIONAL METHODOLOGY 
A. Physical quantities (selfen_phon, selfen_elec, a2f) 

In this section we describe some of the physical quantities which can be calculated using EPW. The imaginary part 
of the phonon self-energy within the Migdal approximation [36l |37| is calculated as: 

q^ =Im 2^ Wk|5„m(k,q)| — — ■ (1) 

^^ Enk - Emk+q " '^<W + «^ 

In Eq. ([T} the electron-phonon matrix elements are given by 

5mn(k,q) = (V'mk+q|9qj.y|-0„k), (2) 

with -^nk the electronic wavefunction for band m, wavevector k, and eigenvalue e„k, da^nV the derivative of the 
self-consistent potential associated with a phonon of wavevector q, branch index v^ and frequency ojqjy. The factors 
/(frik), /(cmk+q) hi Eq. ([T|) are the Fermi occupations, and Wk are the weights of the k-points normalized to 2 in 
order to account for the spin degeneracy in spin-unpolarized calculations. A very common approximation to Eq. ([1]) 
consists of neglecting the phonon frequencies Lo^y and taking the limit of small broadening r\. The final expression 
is positive definite and is often referred to as the "double-delta function" approximation [36|. This approximation 
is no longer necessary when using EPW. We note that the imaginary part of the phonon self-energy in Eq. ([1]) also 
corresponds to the phonon half-width at half- maximum 7qy. 

The electron-phonon coupling strength associated with a specific phonon mode and wavevector Aq^ is given by 

Aqi. = — V w;k|g™„(k,q)p(5(e„k)(5(e„k+q), (3) 

JVp^qi/ — ' 

with 5 being the Dirac delta function. In the double-delta function approximation the coupling strength Xqi, can be 
related to the imaginary part of the phonon self-energy 11" as follows: 

1 n" 

V-^^ (4) 

The total electron-phonon coupling A is calculated as the Brillouin-zone average of the mode-resolved coupling 
strengths Aqj^: 



A^^WqAq,.. (5) 



qi/ 



In Eq. ([5]) the Wq are the Brillouin zone weights associated with the phonon wavevectors q, normalized to 1 in 
the Brillouin zone. The Eliashberg spectral function a^F can be calculated in terms of the mode- resolved coupling 
strengths Aq^ and the phonon frequencies using: 

a^F(a;) = -^u;qa;q^Aq^(5(w -cjq^). (6) 

qiy 



The transport spectral function o^Ft is obtained from the EUashberg spectral function a^F by replacing Aq^ with 

AT,qi. = -r; V Wk|gj;i„(k,q)p(5(e„k)<5(emk+q) I 1 —, ^^^^ I , (8) 

Nptdqi, ^ V v„kr / 

with v„k = VkEnk the electron velocity. 

The real and imaginary parts of the electron self-energy I]„k — Sj^j, + JS"^. can be calculated as 



Snk = 2Xl^ql5mn(k,q) 
qv 

with n(iOq^) the Bose occupation factors. 



njuJci;,) + /(£mk+q) njuq,,) + 1 - /(e„tk+q) 



Cjik — Cmk+q + l^qi/ " iV Ejik — fmk+q — ^qf + «'7 



(9) 



B. Calculation of the matrix elements on the coarse Brillouin zone grid (elphon_shuf f le_wrap) 

The key task of EPW is to calculate from first-principles electron-phonon matrix elements for a large number of 
electron and phonon wavevectors with a modest computational effort. The initial step of this process is to determine 
the electron-phonon matrix elements on coarse grids of electron and phonon wavevectors using standard DFT and 
DFPT methods 38J. The result of this step is a set of matrix elements given by Eq. ^. The wavefunctions and the 
phonon perturbation potentials are read into EPW, then the matrix elements on the coarse grid are calculated within 
EPW. 

The phonon dynamical matrices and the linear variations of the self-consistent potential are read from file only 
for wavevectors in the irreducible wedge of the Brillouin zone. The dynamical matrices and the potentials variations 
for all the remaining points of the uniform grid are generated using crystal symmetry operations. This strategy is 
advantageous since the computation of the phonon dynamical matrices and potential variations is generally the most 
time-consuming part of an electron-phonon calculation. For example, a system with cubic symmetry requiring a 
coarse mesh of 8 x 8 x 8 phonon wavevectors needs only have phonons calculated at 29 irreducible points; hence the 
reduction in computational time is a factor of 512/29. 

The electron-phonon matrix elements and the dynamical matrices thus calculated at each point of the coarse 
Brillouin-zone grid can be written to disk through an input file option (prefix. epb files). Subsequent executions of 
EPW can forgo the recalculation of the matrix elements and dynamical matrices by reading these data from disk. 

C. Interface with wannier90 (pw2wan90epw) 

wannier90 generates maximally localized Wannier functions by minimizing the spread of the Berry-phase posi- 
tion operator through a unitary transform [16|. Details on how to run wannier90 can be found in the wannierQO 
documentation [2i]. 

In standalone mode wannierQO reads three files generated from a first-principles calculation (prefix. mmn, 
prefix. cunn, prefix. eig), a file which defines the starting guess for determining MLWFs and the crystal struc- 
ture (pref ix.nnkp), and a runtime input file (prefix. win). The execution of WcinnierQO therefore requires the user 
to run multiple programs and handle the files to be passed between wannier90 and Quantum-ESPRESSO. In order 
to simplify the calculation of the electron-phonon coupling EPW calls WcinnierQO as a library. EPW therefore requires 
only the wavefunction files from Quantum-ESPRESSO and a runtime input file in order to determine MLWFs using 
wannier90. The quantites required for running wannierQO are either calculated within EPW (for instance the overlap 
matrices A^^ and M,„„ [2]), or else read in from file (for instance the Kohn-Sham eigenvalues). These data are 
then passed from EPW to wannier90 through the wannierjrun library routine. This feature of EPW ensures that the 
execution of wannier90 is embedded within EPW. Hence the end user is only required to run one single executable 
which communicates directly with wannier90. EPW also includes an option to pass additional input parameters to 
wannier90. This allows the user to access all the available features of wannier90, such as for instance plotting 



bandstructures or MLWFs. When called from EPW, wannier90 produces the gauge matrix [/^„ [2|, which yields the 
transformation between Bloch eigenstates and MWLFs, according to: 



n 



(27r) 



WnR^ (r) = -—r dk e-*''-^= > ^ [/^„ Vmk(r), (10) 



BZ 



/ J van 



where WnR^ is a MLWF associated with the direct lattice vector Re, i7 is the unit cell volume, and the integral is 
discretized over the Brillouin zone. The array [/^„ has the dimensions of the number of Bloch bands times the number 
of MLWFs times the number of electronic wavevectors in the coarse Brillouin zone grid. This matrix is written to 
disk and can directly be read on subsequent program executions. The array C/^„ is used in EPW in order to transform 
Bloch functions into MLWFs. 

D. Transformation from the Bloch representation on the coarse Brillouin zone grid to the MLWF 

representation (ephwaim_shuff le) 

After calculating the electron-phonon matrix elements in the Bloch representation for each wavevector on the coarse 
electron and phonon Brillouin zone grids, EPW transforms the electronic Hamiltonian, the phonon dynamical matrix, 
and the electron-phonon matrix elements into the Wannier representation. Derivations and detailed explanations of 



the following transformations can be found in Ref. jl5{ . For clarity the electron band and the phonon mode indices will 
be omitted in the following equations. The electronic Hamiltonian in the MLWF representation H^ ^, is obtained 



as: 



<,R^, - E^ke-''(^-^^)L/^i/£'[/k. (11) 

k 

In this case the weights w\^ are normalized to 1. The Hamiltonian matrix elements in the Wannier representation 
H^ R' decay rapidly with the distance |Ro — R^l, as they scale with the overlap of MLWFs centered in the unit cells 
R = Re and R = RJ,, respectively. The dynamical matrix can be transformed to a localized real-space representation 
using 

<,H^ = E «^qe-^'»(^;-^p)e,<e^, (12) 

q 

where Sq are the orthonormal eigenvectors of the dynamical matrix. The matrix D^ j^, is the mass-scaled matrix of 
force constants. The electron-phonon matrix elements in the MLWF representation are obtained using: 

ff(Re,Rp) = -^5]«;k«^qe-'(''-^=+'i-^p);7^+q5(k,q)f/kUqi (13) 



k,q 



where the Uq^ are the phonon eigenvectors scaled by the atomic masses [15J. In order to check the spatial decay of 

H^ Q, D^ Q, and (7(Rc,Rp) the magnitude of these quantites as a function of Re and Rp is written to formatted 
files in the working directory. A run-time option allows for all data in the Wannier representation to be written into 
one single file. For subsequent calculations these data can be read in and program execution can restart without the 
need to go through the prior computational steps. 

E. Tranformation from the MLWF representation to the Bloch representation on the fine Brilluin zone grid 

(ephwaim_shuf f le ) 

The accuracy of EPW calculations depends on the spatial localization of the MLWFs and the phonon perturbations. 
Typically MLWFs are localized within a few A[39|, |40|. The localization of the phonon perturbation is dependent 
upon the dielectric screening properties of the system, and must be verified in each case before proceeding with the 
interpolation [l5[. 

In order to calculate the electronic eigenstates, the phonon frequencies, and the electron-phonon matrix elements 
on a fine Brillouin zone grid, the Hamiltonian, the dynamical matrix, and the electron-phonon matrix elements are 
truncated outside of a real-space supercell containing A^k and A^q unit cells in the case of electrons and phonons, 



respectively. Here A'k and A'q are the number of grid points in the coarse BriUouin zone meshes. Following this 
truncation it is possible to perform an interpolation back into the Bloch representation onto arbitrary electron and 
phonon wavevectors. The Hamiltonian, the dynamical matrix, and the electron-phonon matrix elements are Fourier- 
transformed back to the Bloch representation by inverting Eqs. (fTTj) . P^ . and ([T^ |l5| . 

The procedure described here enables the calculation of electron-phonon matrix elements on extremely fine BriUouin 
zone grids. The fine grids of electron and phonon wavevectors are specified by the user in input. The procedure 



adopted here is similar to the Fourier interpolation of time series commonly adopted in signal processing |41| . The 
same strategy as the one employed in EPW has been adopted in order to study Fermi surfaces tlq and the anomalous 
Hall effect 0. 

At this stage the Hamiltonian, the dynamical matrix, and the electron-phonon matrix elements on fine grids of 
electron and phonon wavevectors are used to compute the physical quantities described in Sec. IHI Al 

F. Summary of the tasks executed by EPW 

The computational steps described in the previous sections can be summarized schematically as follows: 

1. Electron eigenstates and eigenvalues are read from disk, 

2. wannier90 input data including the overlap matrices A,„„ and Mmn are computed, 

3. wannier90 is called as a library, and the resulting MLWF localization matrix is stored on disk. 

4. The electron-phonon matrix elements are computed on a coarse grids of electron and phonon wavevectors in the 

Bloch representation, 

5. The Hamiltonian, dynamical matrix, and electron-phonon matrix elements are transformed from the Bloch repre- 

sentation to the MLWF representation, 

6. The Hamiltonian, dynamical matrix, and electron-phonon matrix elements are Fourier-transformed back to the 

Bloch representation on arbitrarily dense grids of electron and phonon wavevectors, 

7. The electronic eigenvalues, phonon frequencies, and electron-phonon matrix elements are processed in order to cal- 

culate physical quantites such as, for instance, the electron self-energy or the electron-phonon coupling strength. 

IV. PARALLELIZATION 

EPW is designed and optimized to be executed on multiple processors. Within the Quajitum-ESPRESSD package, 
parallel tasks are split into processor "pools" , where each pool is allocated a set of electronic wavevectors k. If there 
is more than one processor within each pool, the reciprocal space G-vectors for describing the wavefunctions are split 
amongst the processors of that pool. EPW can be executed in parallel by splitting the calculation over electron or 
phonon wavevectors. EPW is not parallelized over the reciprocal-space G-vectors. The parallclization strategy in EPW 
is tailored to each step of the computation. As execution of EPW begins, the coarse k-point mesh is generated and 
distributed amongst the available pools, as in Quantum-ESPRESSD. 

Most of the inputs to the wannier90 such as the lattice and reciprocal vectors are passed through Quantum-ESPRESSD 
and EPW. Important exceptions are the overlap matrices Amn and Mmn- Within EPW the computation of the overlap 
matrix elements is distributed across the coarse grid k-point pools. As the elements of the matrices Amn and Mmn 
are computed independently, the execution time of this step scales inversely with the the number of pools. 

wannier90 is executed serially on each processor through a library call. The time required to determine MLWFs is 
negligible as compared to the computation of the electron-phonon matrix elements on the coarse BriUouin zone grids. 

The electron-phonon matrix elements on the coarse BriUouin zone grid are computed sequentially for each phonon 
wavevector, while the electron wavevectors are distributed across pools. 

The computation of the electron-phonon matrix elements on the fine BriUouin zone grids is parallelized over elec- 
tronic or phonon wavevectors depending on the calculation type. For example, when calculating the phonon self-energy 
an integration is required over electronic wavevectors. In this case the most convenient strategy is to parallelize over 
the electronic wavevectors. The converse is true for the calculation of the electron self-energy. In those cases where 
the final quantity requires integrations over both electron and phonon wavevectors, such as the total electron-phonon 
coupling strength, the efficiency of EPW is independent of the choice of the parallclization scheme. 
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FIG. 1: Parallelization of EPW: speedup vs number of processors observed in parallelizing a) the generation of MLWFs, b) the 
calculation of the electron-phonon matrix elements on the coarse Brillouin zone grids, c) the interpolation of the electron-phonon 
matrix elements to the fine Brillouin zone grids as well as the computation of the phonon self-energy. Panel d) displays the 
speedup observed to perform the entire computation including the initialization step and the intermediate I/O. The speedup is 
defined as the ratio between the time it takes to run a calculation on a single processor and the time it takes to run the same 
calculation on a given number of processors. The diagonal dotted lines correspond to the ideal speedup, which is equal to the 
number of processors employed. The calculations have been performed for hole-doped SiC in the rigid-band approximation. 
The unit cell contains 2 atoms. The blue line corresponds to a calculation with coarse Brillouin zone grids containing 6x6x6 
k- and q-points, a kinetic energy cutoff of 30 Ry, and fine Brillouin zone grids containing 1000 points each. The red line 
corresponds to a calculation with coarse Brillouin zone grids containing 8 x 8 x 8 k- and q-points, a kinetic energy cutoff of 60 
Ry, and fine k- and q-points grids containing lO'' and 104 points (the irreducible points of a uniform and unshifted 14 x 14 x 14 
grid), respectively. The parallelization becomes more efficient as the number of electron and phonon wavevectors in the fine 
Brillouin zone grids is increased. 

Figure [T] shows the efficiency of our parallelization strategy for various sections of EPW for the test case of hole-doped 
SiC. This example is included for reference in the EPW distribution. 

A. Calculations "on the fly" 

The interpolation part of EPW can be executed in two different modes. In the first mode EPW calculates all the 
electron-phonon matrix elements for every k- and q-point of the fine Brillouin zone grids, and stores this information 
to disk. The subroutines for calculating the electron or phonon self-energy are then called and the matrix elements 
along with the electronic eigenvalues and the phonon frequencies are read from disk in order to evaluate Eqs. ([T]), ([9]). 

In the second mode of operation, the electron or the phonon self-energy are calculated "on the fly" as the matrix 
elements are generated, and all the quantities are overwritten. In this case the matrix elements are not stored on disk. 

This feature allows us to address systems for which the storage required for the electron-phonon matrix elements 
would be exceedingly large. For example, if the fine electron and phonon Brillouin zone grids consisted each of 10^ 
wavevectors, and we had 8 MLWFs and 6 phonon modes, then the entire set of double precision matrix elements 
would require more than 1 TB of storage. On the other hand the calculation on the fly would reduce the storage 
requirement down to a manageable size of 1 TB. 

V. EXAMPLES 

In the following we illustrate the capabilities of EPW by describing three prototypical systems: (i) Lead, a simple 
metal which is also a strong-coupling superconductor, (ii) Graphene, a two-dimensional zero-gap semiconductor with 
linear quasiparticle dispersions close to the Fermi level, (iii) Silicon carbide, a wide-gap semiconductor which becomes 
a low-temperature superconductor upon boron doping. 



Lead 



Solid Pb is the prototypical strong-coupling superconductor. While the superconducting transition temperature of 
Pb is less than 10 K, the electron-phonon coupling strength has been measured and calculated to be very large, in 
the range of A = 1.3 ± 0.1 ^^. 
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FIG. 2: Spatial decay of the largest components of the Hamiltonian H^^ q (a), the dynamical matrix D^ q (b), and the 

electron-phonon matrix elements (;(0, Rp) (c) and gr(Re,0) (d) for fee Pb. The data are plotted as a function of distance and 
are taken along several directions. 



In this section we present an example calculation of the electron-phonon coupling in bulk Pb and we compare 
our results to experimental data. The calculations were performed within the local-density approximation to density 
functional theory using Quantum-ESPRESSO, wannier90, and EPW. A norm-conserving, scalar-relativistic pseudopoten- 
tial including four valence electrons was used to take into account the core-valence interaction. The electronic states 
were computed within a plane-wave basis with a kinetic energy cutoff of 80 Ry. The charge density was computed 
self-consisently on a 16 x 16 x 16 F-centered Brillouin zone grid. The electronic wavefunctions used within EPW were 
calculated on an 8 x 8 x 8 uniform Brillouin zone grid. The phonon dynamical matrices and the linear variations 
of the self-consistent potential were calculated within DFPT using the Quantum-ESPRESSO package, using the same 
convergence parameters as above. We considered a uniform 8x8x8 Brillouin zone grid for the phonon calculations, 
corresponding to 29 irreducible phonon wavevectors. 

Four Wannier states were used to reconstruct the electronic structure near the Fermi level. These states were found 
to be sp'^-like functions localized symmetrically along each of the Pb-Pb bonds, with a spatial spread of 2.17 A. The 
spatial decay of the electron Hamiltonian, the phonon dynamical matrix, and the electron-phonon matrix elements in 
the MLWF representation are shown in Fig. [2l 

In order to calculate the phonon linewidths and the total electron-phonon coupling, the electron Hamiltonian, the 
phonon dynamical matrix, and the electron-phonon matrix elements were transformed from the MLWF representation 
to the Bloch representation by the generalized Fourier interpolation described in Sec. lIIIEl By using 10^ k-points, 8000 
q-points, and a smearing parameter of 30 meV, we obtained a total electron-phonon coupling strength A = 1.41. The 
Eliashberg spectral function a'^F{io) is shown in Fig. |3l together with the experimental curve obtained by inverting 
tunneling data [42| . Our method clearly yields a very good agreement with experiment. 

For comparison with experiment we also calculated the superconducting transition temperature Tc starting from 
the total electron-phonon coupling. We adopted the modified McMillan equation J47| and a Coulomb repulsion 
parameter /i* between 0.10 and 0.15, obtaining Tc equal to 4K to 6K respectively. This result compares favorably to 
the experimental value 7.2 K. 
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FIG. 3: Eliashberg spectral function of fee Pb. The solid line is from Ref. [42], the dashed line has been obtained using EPW. 



B. Graphene 



Graphene is a two-dimensional sheet of carbon atoms which has attracted enourmous attention from the scientific 
community due to its peculiar electronic properties [48| , and in particular the Dirac-like nature of electrons near the 
Fermi level 49|, |50[. Considerable attention has been paid to the measurement of many-body renormalization effects, 
such as the effect of the electron-phonon interaction on the electronic bandstructure, using for instance angle-resolved 
photoemission experiments [51] • 

We here present the results of an electron-phonon calculation on a sheet of freestanding graphene [23| . The local 
density approximation to density functional theory was employed in conjunction with a norm-conserving carbon 
pseudopotential. A plane wave basis with a kinetic energy cutoff of 60 Ry was used. The graphene sheets were 
separated by a vacuum of 8 A in a supercell geometry in order to eliminate spurious interactions between periodic 
replicas. The relaxed C-C bond length was found to be 1.405 A. 

Three in-plane bonding Wannier functions, two Pz-like Wannier functions (one per each carbon atom), and two 
s-like Wannier functions directly above and below the center of the hexagon (away from the graphene plane by 1.57 
A) in the two-atom unit cell were used to describe the electronic structure. The spatial spread of these MLWFs are 
0.565 A, 0.782 A, and 1.726 A, respectively. The spatial decay of the electron Hamiltonian, the phonon dynamical 
matrix, and the electron-phonon matrix elements in the MLWF representation for a 6 x 6x 1 grid are shown in Fig. 

SI 

The electron self-energy was computed and the integral of of Eq. ^ was performed with 10^ phonon q-vectors, 
obtained through the interpolation method of this communication. The mass renormalization parameter, A„k has 
been obtained as an energy derivative of the electron self-energy as in A„k — ^R-^ S„k(enk)|e=eF- The results of 
the calculations of the real and imaginary parts of the electron self-energy as well as the electron-phonon coupling 
parameter as a function of Fermi energy are given in Fig. [Sj 



Silicon carbide 



The possibility of achieving superconducting behavior in semiconductors has recently attracted considerable inter- 
est, and several experimental 52|, [53 1 and theoretical 2^, 127], |5J-l59| studies have been performed on this class of 
materials. Boron-doped SiC is a promising candidate because of its potential uses in power electronics owing to its 
large breakdown voltage. 

For this test case we performed electron-phonon calculations on a rigid-band model of 4% hole-doped cubic SiC. 
For our calculations we employed the local density approximation to density functional theory, norm-conserving 
pseudopotentials, and a plane wave basis with a kinetic energy cutoff of 60 Ry. The relaxed lattice parameter was 
found to be 4.34 A, in good agreement with the experimental value of 4.36 A. EPW was executed using coarse meshes of 
8x8x8 uniform grids of electron and phonon wavevectors. Four Wannier functions per formula unit were considered 
to describe the valence electronic struture. MLWFs were found to be sps-like functions with an average spread of 
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FIG. 4: Spatial decay of the largest components of the Hamiltonian H^^q (a), the dynamical matrix D^ q (b), and the 
electron-phonon matrix elements g{0, Rp) (c) and g(Ro, 0) (d) for graphene. The data are plotted as a function of distance and 
are taken along several directions. 
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FIG. 5: The real (a) and imaginary (b) parts of the electron self-energy Ek near the Dirac point as well as the electron-phonon 
coupling A„k as a function of doping (c) for graphene. In panels (a) and (b), the Fermi energy has been set to zero. 

1.15 A. The spatial decay of the electron Hamiltonian, the phonon dynamical matrix, and the electron-phonon matrix 
elements in the MLWF representation are shown in Fig. [51 We calculated the total electron-phonon coupling strength 
A — 0.34 using 250,000 k-points and 8000 q-points. The superconducting transition temperature was calculated using 
a Coulomb repulsion parameter fi* = 0.1 and was found to be Tc = 1 K. Our result is in good agreement with recent 
experimental data yielding Tc = 1.4 K|52|. 



VI. CONCLUSION 



In this communication we introduced EPW, a computer code for calculating the electron-phonon coupling from 
first-principles using density functional perturbation theory and maximally localized Wannier functions. EPW enables 
extremely accurate and highly efficient calculations of the mode-resolved and total electron-phonon coupling strength, 
as well as phonon and electron linewidths. The code is distributed through the website http://epw.org.uk and is 
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FIG. 6: Spatial decay of the largest components of the Hamiltonian -ffR^_o (&)i the dynamical matrix D^ q (b), and the 
electron-phonon matrix elements (;(0, Rp) (c) and (?(Rc,0) (d) for cubic hole-doped SiC. The data are plotted as a function of 
distance and are taken along several directions. 



available under the terms of the GNU General Public License. Plans are in place to extend EPW in order to implement 
the anisotropic Eliashberg theory [60| , the density functional theory for superconductors 6l|, |62| , and phonon-assisted 
optical responses [63|. 



VII. APPENDIX 

A. Specifying a unique gauge for the electronic wavefunctions (setphases) 

Since nondegenerate electronic wavefunctions are uniquely defined up to a phase, and a set of degenerate wave- 
functions can be mixed via a unitary matrix, electron-phonon matrix elements are machine dependent. In some cases 
it may be convenient to use electron-phonon matrix elements outside of EPW, for instance in the study of phonon 
sidebands in excitonic spectra [6J], or in the phonon-assisted Auger recombination [65|. In order to meaningfully use 
the electron-phonon matrix elements generated by EPW outside of the code it is necessary to define uniquely the phase 
of each wavefunction as well as the way in which degenerate wavefunctions are mixed. 

In order to set the gauge of each wavefunction in EPW we proceed as follows. First, we determine the subset of degen- 
erate wavefunctions at each wavevector on the coarse mesh. Then we artificially lift the degeneracies by diagonalizing 
the subset of degenerate states with respect to a small external perturbation. In principle the perturbation could take 
on any form, however for convenience EPW employs a nonlinear combination of the data found in the prefix. dvscf 
files. Finally, the Fourier components of each wavefunction are scaled by a complex factor exp{i0) in such a way 
that the largest component of each wavefunction is real-valued. In this way, the machine dependence of the matrix 
elements is eliminated. This procedure is explained in detail in Ref. [15|. 
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